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A scalar theory of gravitation with a preferred reference frame (PRF) is considered, that accounts 
for special relativity and reduces to it if the gravitational field cancels. The gravitating system 
consists of a finite number of perfect- fluid bodies. An « asymptotic » post-Newtonian (PN) 
approximation scheme is used, allowing an explicit weak-field limit with all fields expanded. Exact 
mass centers are defined and their exact equations of motion are derived. The PN expansion of 
these equations is obtained: the zero-order equations are those of Newtonian gravity (NG), and the 
equations for the first-order (PN) corrections depend linearly on the PN fields. For PN corrections 
to the motion of the mass centers, especially in the solar system, one may assume « very-well- 
separated » rigidly moving bodies with spherical self-fields of the zero-order approximation. The 
PN corrections reduce then to a time integration and include spin effects, which might be 
significant. It is shown that the Newtonian masses are not correct zero-order masses for the PN 
calculations. An algorithm is proposed, in order to minimize the residual and to assess the velocity 
in the PRF. 

1. Introduction 

In 1905, Einstein suggested that the very notion of an absolute motion does not make sense in the 
physical world, and abolished the ether as « superfluous. » However, the « principle of relativity » 
was formulated in 1904 under this same name by Poincare [1], in the following form: « the laws of 
physical phenomenons must be the same, whether for a fixed observer, as also for one dragged in a 
motion of uniform translation, so that we do not and cannot have any mean to discern whether or 
not we are dragged in a such motion. » Thus, for Poincare, the notion of an absolute motion does 
make sense, even though « it seems that this impossibility to demonstrate the absolute motion be a 
general law in Nature » [2]. Moreover, Poincare, as well as Lorentz, always reserved the possibility 
that the principle of relativity might be falsified by an observation, and indeed both considered an 
ether until the end of their lives (see e.g. Ref. 3). 

The principle of relativity, as formulated either by Poincare or by Einstein, leads to the same 
physical theory of special relativity (SR), which is one of the best-verified theories of physics. But 
SR does not involve gravitation. Although Einstein's theory of gravitation, which is the currently 
accepted theory of gravitation (see e.g. Will [4]), is named « general relativity » (GR), it is well- 
known since Fock [5] that the principle of relativity does not hold true in GR. For instance, it is 
clear that an observer can discern in a number of ways whether or not she is in a motion with 
respect to the static reference frame of Schwarzschild's space-time. It is true that, due to its 



manifest covariance under arbitrary coordinate changes, GR does not have any a priori preferred 
reference frame, in the following sense: if, in GR, we model the solar system (say) as isolated, we 
do not have to worry about the motion of the solar system with respect to whatever reference frame, 
and thus may choose our coordinate system such that, for instance, the Sun is at rest in that 
coordinate system. (In reality, however, the solar system is not isolated: it should be embedded in a 
large-scale model, thus the metric in the solar system should be « matched » with, say, the 
Robertson- Walker space-time metric representing the assumed macroscopically-homogeneous 
Universe.) But the point is that the experimental support of SR does not, as such, provide much 
constraint on the theory of gravitation: one may define a very large class of « relativistic » 
modifications of Newton's gravity (NG), by the condition that any of them should reduce to SR 
when the gravitational field cancels. (In GR, following Synge [6], this may be defined as the case 
where the Riemann tensor cancels, in which case the Minkowski space-time is indeed obtained.) If 
one defines in this way the « relativistic » character of a theory of gravitation, nothing forbids that a 
such theory may have an a priori privileged reference frame - thereby violating, not only the 
explicit form of the principle of relativity (as is already the case for GR), but its very spirit. If a 
privileged reference frame would manifest itself in the presence of a gravitational field (while 
remaining « hidden » in its absence), it would be easier to extend quantum mechanics to the 
situation with gravitation [7] . 

Many consider that, among relativistic theories of gravitation, it is unlikely that any 
preferred-frame theory could make sensible predictions. However, it has recently been proposed [8- 
10] a theory which, despite the simplicity which is conveyed to it by the fact that it is a scalar 
theory, and despite its preferred-frame character, has the following encouraging four properties: (i) 
it is relativistic, i.e., it does reduce to SR when the gravitational field cancels, (ii) It has the correct 
Newtonian limit [11]. More precisely, one may define, within that theory, a small parameter £ that 
characterizes the weakness of the gravitational field in the considered system and that is 
approximately equal to the ratio of the maximum orbital velocity in the system, to the velocity of 
light. And one may build an asymptotic scheme in which, up to the order e 2 not included, the 
expanded equations reduce to NG [12-13]. This means that this theory is (at least) as good as NG in 
the solar system, (iii) It predicts [15] the observed acceleration of the cosmic expansion [16, 17] 
(iv) It recovers exactly the standard post-Newtonian (PN) predictions of GR as to the gravitational 
effects on light rays [14]. This is a crucial point, because those effects represent the best- verified 
consequences of GR, and nowadays the most important ones practically. Yet it is not enough that 
this theory be as good as NG for celestial mechanics: one expects from a viable theory of 
gravitation that it improve celestial mechanics as compared with NG. This requirement is usually 
concentrated on the demand that Mercury's residual advance in perihelion should be explained. 
However, the modern process of observational test of the celestial mechanics based on GR involves 
a general least-squares fitting of a large set of observational data by using the equations of that 
celestial mechanics [18]. In accordance with this process, the assessment of different celestial 
mechanics should not be based on God's judgment: « are Mercury's 43" predicted or not ? » but, 
more prosaically (and more accurately), on the comparison between the respective least-squares 
residuals. In doing so, one may make different comparisons based on different choices of the set of 
input data. E.g., one may check how well-predicted is Mercury's perihelion motion, by not 
incorporating it in the input data. 
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The aim of this paper is to make an important theoretical step in order to make possible the 
numerical implementation of a celestial mechanics based on the scalar theory in an astronomical 
ephemeris. This step consists in obtaining the equations for the first PN corrections to the motion of 
the mass centers and in putting them in a tractable form. We shall base this work on foregoing 
results [13] that give the local equations of motion and the boundary conditions in the first PN 
approximation (PNA), hence Section 2 summarizes these results. We emphasize that an 
« asymptotic » PNA is considered, in the sense that the usual method of asymptotic expansion is 
applied to a one-parameter family of similar gravitating systems, constructed from the data of the 
physical system of interest [13]; as a consequence, the gravitational field and the matter fields are 
expanded, as in Refs. 19 and 20 - in contrast with the standard PNA [4-5, 21-23], in which only the 
gravitational field is expanded. Moreover, we consider extended bodies, as is relevant to celestial 
mechanics [24]. To our knowledge, this is the first time that an « asymptotic » PNA is used for 
extended bodies. According to Asada & Futamase [25], «the theory of extended object in general 
relativity is still in preliminary stage for the application to realistic systems. » In Section 3, we 
describe how to properly define mass centers in order that their motion be relevant to the 
comparison with astronomical observations. In Section 4, we give the general form of the equations 
for the PN corrections to the motion of the mass centers. In Section 5, we consider the case of well- 
separated bodies, that introduces the maximum ratio of the size of the bodies to the distance 
between bodies as a small parameter 77. In Section 6, we consider the case of « very-well- 
separated » bodies, in which case only terms of the lowest order in 77 are conserved, together with 
the further relevant simplification that occurs when one assumes spherically symmetric zero-order 
self-fields (only when evaluating the PN corrections). In Section 7, we propose an algorithm for the 
numerical implementation of the method and we show that the Newtonian masses must be 
corrected. Some integrals are computed in Appendices. 



2. Resume of previous results 

The theory of gravity which is considered is a scalar theory with a preferred frame E (« ether ») [8], 
with two space-time metrics: a flat « background metric » y° and a curved « physical metric » y, the 
relation between the two being defined through the scalar field [9], and with the motion being 
governed by an extension of Newton's second law rather than by Einstein's geodesic assumption 
[10]. A brief summary of the theory is given in Ref. 13 (a more detailed one is given in Ref. 14; 
Ref. 26 presents the complete construction of the theory, though without detailed calculations). 

2.1 Equations for the exact fields 
The preferred frame E is assumed to be an inertial frame for the flat metric, thus there are Galilean 
coordinates (x^) for y° [i.e., (/^v) = (f]nv) = diag(l, -1, -1, -1)], that are bound to the frame E. The 
inertial time T = x°/c in the preferred frame is called the « absolute time ». In any coordinates (y^) 
bound to E and such that y Q = cT, the space-time metric is given as function of the scalar field/ by: 

700=/, Ya = -9ip 7or = 0, (2.1) 

where g is the spatial part of the « physical » metric y in the frame E, which is related to the the 
spatial part g° of the flat metric y° in the frame E as follows: 
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(2.2) 



In the case that/= 1 in some region of space-time (which means that no gravitational field is there), 
the physical metric coincides with the flat metric in that region. Except for cosmological problems, 
the field equation [8-9] may be written as [10]: 



A/- 



/ 



8ti G 



(f = cT), 



(2.3) 



with A the usual Laplace operator, defined with the Euclidean metric g , and where G is Newton's 
gravitation constant, and o = (T ) E is the material mass-energy density in the ether frame. This 
equation is merely space-covariant. Here, we shall consider that each body is made of a perfect 
fluid; that is, the mass-energy tensor of matter and non-gravitational fields has the form 



T» v = (p* + p/c 2 )U»U V - (p/c 2 ) f 



(2.4) 



where t/ M = dx^/ds is the four- velocity of the fluid, p is the pressure and (}^ v ) is the inverse matrix 
of (Yfiv), and where n* is the volume density of the rest-mass plus internal energy in the proper 
frame, expressed in mass units, 

/j*sp*(l +nic 2 ). (2.5) 
The dynamical equations are then given in Cartesian coordinates (such that g°ij= dy) by [13]: 

d T (yu i ) + dj(\i/u i u j ) + r i j k yu j u k + t i k \i/u k = \i/fg i -g ij P j , (2.6) 



d T (xiff-(p/c 2 )) + dj(xi/fu j ) = (o/2)d T f. 



(2.7) 



where 



dx* 
dT 



tf j = (g-%, V jk = ^(gy, k + g lkJ - g jK 0, t l k = \g ij d T g jk , y/^a + p/ (c 2 f). (2.8) 



The material mass-energy density is given by 



0=(T ) E = 



f 



c L J 



7, 



(2.9) 



where y v = (1— v Ic )~ is the Lorentz factor, with v = g^ v v 1 and V = dx l /dt x = uNf. One may 
check that, if the scalar field /can be considered equal to 1 (no gravitational field), then Eqs. (2.6) 
and (2.7) reduce to the dynamical equations for a perfect fluid in SR [5, 22]. However, due to Eq. 
(2.3), the gravitational field can hardly cancel unless there is no matter - as in GR. 
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2.2 Weak-gravitational-field limit 
This is defined in the following way [13]: we consider the gravitating system of interest, S, made of 
perfect fluids and assumed to obey the equations of the scalar ether-theory, Eqs. (2.3), (2.6) and 
(2.7). We define 

V=c 2 (l-/)/2, (2.10) 
V max = max { V(x, 7=0); xeM}, £ = (Wc 2 ) 1/2 = [(l-/ min )/2] 1/2 (2.11) 

(here M is the space manifold, made of all points that are bound to the ether frame). We assume 
that £o is well-defined, finite, and small: this is the first physical assumption about system S, it 
means that the gravitational field is weak. We define a one-parameter family (S e ) of gravitating 
systems by the fields/ 6 , p e , u e , those being defined as the solution (assumed unique for all xe M and 
in some interval [0,72(£)] for 7) of the initial-value problem 

, ? V(x,0) df £ 2E 3 dV 
/ £ (x,0) = l-2<^ 2 ^fA (x,0) = -^V — (x,0), (2.12) 

c oT c dT 

p £ (x, 0) = ?p(x, 0), u e (x, 0) = £u(x, 0), (2.13) 

for equations (2.3), (2.6) and (2.7). Here £ = £/£q and the state equation in system S e is 

P* e = F e (p e ) = ?F(%- 4 p) [or F e (p) ee e 2 F Y (e - 4 p)] , (2.14) 

where p* = F(p) is the (barotropic) state equation in system S. This definition by an initial- value 
problem is justified by the fact that the natural boundary-value problem in the scalar ether-theory is 
indeed the full initial-value problem [13]. 

2.3 Asymptotic expansions of fields and equations 
One expects that the family of fields just defined admits an asymptotic expansion as the small 
parameter e tends towards zero. More precisely, if one changes units for system S e , multiplying the 
starting time unit by £ 1 and the mass unit by e 2 , then £ becomes proportional to 1/c and all matter 
fields are of order zero with respect to e. This is equivalent to say that, in invariable units, one has 

/=ord(e 4 ), u e =ord(e), p* e =ord(£ 2 ), (2.15) 

as in the weak-field limit of NG, and as is implied by the initial conditions (2.13) with (2.14). Since 

2 2 

it is only 1/c that enters the field equations (2.3), (2.6) and (2.7), one states expansions in 1/c for 
the independent fields, adopting the varying units just described: 

f £ = 1-2 U/c 2 - 2 A/c 4 + 0(£ 6 ), (2.16) 
p e = p + pi/c 2 + O(£ 4 ), (2.17) 
u e =u +u 1 /c 2 + O(e 4 ). (2.18) 

[Theoretically, these expansions should be mathematically deduced from the definition (2.12)- 
(2.13).] From these expansions, one then deduces expansions for the other fields, and for the 
equations. Those « expanded equations » are valid in any units, provided (of course) one expresses 



5 



the expanded fields in these units. As to the expansions, they take a slightly different form in 
invariable units [13]. However, we shall consider mainly the 1PN approximations of the fields, thus 
omitting the remainder terms in the expansion: e.g. 

p e {l) = po + Pl / c 2 . (2.19) 

Those keep the same expression in any units although, in invariable units, po is like £ 4 and p\ like 
£ 6 . Thus we shall use fixed units in this paper, unless otherwise mentioned. The order in £ of an 
expression is easily obtained by passing through the varying units (in which the order is directly 
read out), but we shall leave it to the reader. Anyhow, we emphasize that, in this work, we simply 
compute all expressions consistently at the first PNA, i.e. the second approximation in the 
expansion in powers of e 2 . We also shall omit the superscript £ for the 1PN approximations of the 
fields. Note that we are actually interested in the given system S, that corresponds to the given, 
small value £o for £. Finally, we shall omit the index for the zero-order coefficients in the 
expansions: thus Eq. (2.19) is rewritten as 

P(i) =P + P\/c 2 . (2.20) 

To avoid confusions, the exact fields, when needed, will now be denoted by the subscript « exact ». 
For instance, the density of rest mass in the preferred frame and with respect to the Euclidean 

12 3 

volume measure dV= dx dx dx (in Cartesian coordinates) is defined by 

Pexact = dmo/d V = P*exact Mexact/vfcact • (2.21) 

We shall need the following relations between the expanded fields [13]: 

p*=F(p), p*i=F'(p).p u (2.22) 
p = p*, pi=p* 1 + p(u 2 /2+^, (2.23) 
y/=c7=p=p*, (7i = pi+p(u 2 /2 + 77+ U), y/i = Oi+p. (2.24) 
The expansion of the equations is as follows [13]: Eq. (2.3) is expanded to 

AU = - 4nGp, (2.25) 
AA = - 47lGc7i + d 2 U/dT 2 ; (2.26) 

the energy equation (2.7) is expanded to: 

d T p + dj(pu>) = 0, (2.27) 
d T px + dj(pxu j + pm j ) = 0; (2.28) 
and the local equation of motion (2.6) is expanded to: 

d T ( pu^ + dji pu^j) = pU { - p j, (2.29) 
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d jlpui +y/ l u l ) + d j (y/ l u l u j + pu 1 'u j + pu l u l j ) + ^ pu J u k + pu j d T k ij = 

= - p u i + 2 kij Pj + yf! U i + p A ,■ - 2C/p[7 , (2.30) 

where 

= ^y,* + kikj ~ k jk i , = Uh tj, Hij = ' ' • (2.31) 

u ,k u ,k 



2.4 Boundary conditions for the expanded fields 
Entering the expansions of the fields into the initial conditions (2.12)-(2.13), one gets initial 
conditions for the PN expansions of the matter fields [13]: 

p(x,T = 0) = pvUxJ = 0), Pl (x,0) = 0, (2.32) 

p(x,0) = F[p exact (x,0)], Pi (x,0) = [p(u 2 /2 + U)](x,0), (2.33) 

u(x,0) = u exa ct(x,0), ui(x,0) = 0. (2.34) 

Since the first- approximation equations: (2.25), (2.27) and (2.29), are the equations of NG, the 
gravitational potential U plays just the role of the Newtonian potential in NG. Hence we assume 
spatial-boundary conditions at infinity ensuring that the solution of (2.25) is indeed the Newtonian 
potential associated with p: 

(Vr ) U(x, T) = 0(l/r) and NU (x, T)\ = 0(l/r 2 ) as r = Ixl -> °°. (2.35) 

This represents a physical condition imposed on the field /, in addition to the requirement that £q 
«1 [13]. In the latter work, we proposed similar conditions for the gravitational potential A, but 
this turns out to be inappropriate. Indeed it is not A but a « part » of it that must behave like a 
Newtonian potential [see Eqs. (4.14) and (4.15) below]. 



3. Definition of the mass centers 



In the « asymptotic » scheme considered here, the decoupling between the Newtonian equations and 
the equations for the PN corrections is particularly clear. Hence we shall have the possibility to 
make stronger assumptions for the calculations of the PN corrections than for the Newtonian 
calculations, because the former are very small. It is in particular possible to assume spherical 
bodies at the level of the PN corrections to the mass centers, while keeping a more general 
description at the Newtonian stage. Moreover, as regards the equations for rotational motion, it is 
possible to content oneself with Newtonian calculations, in which case one merely has to evaluate 
the PN corrections to the motion of the mass centers, indeed. However, whereas in Newtonian 
celestial mechanics, only the distribution of the (Newtonian) mass influences the motion of the 
mass centers, the mass-energy equivalence changes this situation for a « relativistic » theory of 
gravitation: i) it is not obvious to decide which mass-energy density is relevant to the definition of 
the mass centers, ii) once a definition has been selected [e.g. one might think to the « active » 
energy density o, Eqs. (2. 3) and (2.9)], one has to expect that it is not only that density which 
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determines the motion of the selected « mass centers ». This means that the internal 
structure (through its energy distribution) and the internal motion (through the corresponding 
kinetic energy) are now a priori expected to play a role. Finally, the « choice » of the energy density 
should be influenced by observational considerations: one should ask if the energy density selected 
is simply related with the emission of electromagnetic radiation detected by the telescope. It seems 
clear that the electromagnetic emission is very loosely correlated with the gravitational energy, 
because the latter is distributed in the whole space, and that instead it is strongly correlated with the 
rest-mass energy. In other words, it is the presence of « matter » in the usual sense, thus 
characterized by its rest mass, that leads to the emission of e.m. radiation. Since we are working in 
the preferred frame of the scalar theory, the rest-mass distribution is defined by the density (2.21), 
whose PN expansion involves the Newtonian density p and the PN correction p\. 

A theoretical argument also leads to taking the rest-mass density (2.21). Indeed suppose one 
defines the « mass centers » as local barycenters for some density <pin some separated domains u) a , 
(i>b , ... [one for each body a, b, the domains depend on time, of course]. Thus one defines the 
position a of the « ^-center » of body a by 

<D fl a = J ©xdV, <D flS J ©dV. (3.1) 



One would wish that the « <p-mass » <i> a be a constant, and that further one have 

<M = I <pu exact dV. (3.2) 



If a continuum moving with velocity w occupies just the domain (d a (T ) at any T, we have 



f d(D I \ 
<S> a=L ^7+div((pwj 



dV, 



(3.3) 



hence also (applying this to £, ' = cp x 1 ): 



-7— f (PxdV=\ <pwdV+{ 



Jt 



+ 



div(<pw) 



dV. 



(3.4) 



In Eqs.(3.3)-(3.4), we may alter the velocity field w(x,T) continuously inside (0 a (T), provided that 
we have %T[(d a (T =0)] = (i> a (T ), where x(T ) = ^t(X) are integral lines of the vector field w (such 
that dx/dr = w on these lines; X is the position at time T =0). Moreover, in the special case that 
(f>(x,T) = on the boundary d(O a (T), an application of the divergence theorem shows that (3.3)-(3.4) 
hold true for an arbitrary vector field w(x,T). Hence, Eqs.(3.1)-(3.2) are certainly satisfied if q> 
obeys the usual continuity equation with the exact velocity u eX act, and there is hope that one may 
substitute some approximate velocity w for u eX act, provided <p and w satisfy the usual continuity 
equation. In the scalar theory, there is an exact conservation equation for the total energy, including 
the gravitational energy [10]. But this conservation equation [obtained by rewriting the r.h.s. of 
(2.7) using the field equation (2.3)] is not the usual continuity equation since it also involves 
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gravitational energy density and flux. Only the rest-mass density obeys the usual continuity 
equation, though it obeys this equation approximately [see Eq.(3. 1 1) below]. 

Thus, we define the exact mass center through the rest-mass density p eX act : 

^aexact = J" fflj PexactXdV, M ?" S h. P^^V . (3.5) 

At the (first) PNA, p exa ct is approximated as p exact = p (V) [1+ 0(e 4 )], with [cf. (2.22), (2.23)] 

p ( i) = p + pi/c 2 , p l = (yx 1 /l+U)F(p)+PiF\p) (3.6) 
Thus, the mass and the mass center are approximated by 

M^ = M a + M l Jc\ M a =\ pdV,M l a = \ Pl dV, (3.7) 

M a ] a (D = L a P (D X dV = M a 3 + M >1 1 

with 

Af.a = J pxdV, Afia,=f p,xdV. 

Note that M a and a are the Newtonian mass and center of mass. Now let us show that M a and M a l 
are (nearly) constant. We first have to precise the definition of the domain (0 a occupied by body a 
(this does not appear to be done in the literature on relativistic celestial mechanics). In the real 
world, there is no sharp separation between different celestial bodies, for a perfect vacuum does not 
exist. This is compatible with different barotropic state equations for different bodies, provided the 
state equations coincide at very low pressure. We define the domain (0 a (T = 0) occupied by body a 
at the initial time by a threshold pressure, thus q(X) = p(X,T = 0) > p m i n a if X belongs to (O a (T = 0). 
Then we follow (O a with the PN velocity as the velocity w : 

w = U(i) = u + u/c 2 . (3.10) 

Now we have by (2.27), (2.28) and (2.15): 

^+div(p (1) u (1) ) = -^div(p 1 u 1 ) = 0(e 7 ), (3.11) 

from where it follows by (3.3) and (3.7)i that M a (l) is conserved at the first PNA: 

dM a (1) /dr=0(e 7 ), (dM a {l) /dT)/M a {1) = 0(e 5 ). (3.12) 

(Actually, a global mass like this is likely to be conserved with an extremely good approximation in 
the scalar theory [27].) Moreover, from the continuity equation (2.27), we get using Eqs.(3.3), (3.7) 
and (3. 10), then (2.15): 

dM a /dT=\ d(SSo P ui .n/c 2 dS = 0(e 5 ). (3.13) 
9 



(3.8) 
(3.9) 



Note that p follows the Newtonian motion that has velocity u ? w. For instance, in the case of 
Mercury, the « Newtonian » body will be, after one century, at 43" angular distance on Mercury's 
orbit from the « PN » body (if the PNA of the theory is good enough to reproduce the observed 
residual advance). This corresponds to some 12.10 km, to be compared with the 2437 km of 
Mercury's « official » equatorial radius. However, by taking /? m in Mercury small enough, the radius of 
Mercury shall be several 10 4 km, say, and yet this « enlarged Mercury » will remain at a 
comfortable distance from other celestial bodies. With this definition, the « Newtonian » position of 
the body and the « PN » one, though at the same distance from one another (roughly 10 4 km after 
one century, in the case of Mercury), will both contain virtually all the matter of the real body, up to 
a negligible amount. In other words, by taking p min a small enough, we may assume that p and pi are 
negligible on d(O a so that, in particular, the rate given by (3.13) [which rate is 0(e 5 ) independently 
of this « enlargement »] is utterly negligible. By (3.12) and (3.7)i, dM a Y /dT is 0(e 5 ) also but, 
actually, from (2.27), we get using Eqs.(3.3), (3.7) 3 and (3.10): 

dM a l /dT =L (pi ui/c 2 - p uO.n dS, (3.14) 

which is utterly negligible with our definition of « enlarged » bodies. But, after some time (may be 
of the order of a century), one should reinitialize » the problem. We get from (3.7)-(3.9): 



(5a = a(i) a = 2 



c 



M a 



-a) 



+ 



o(e 4 ). (3.15) 



As to Ma, we obtain by using (2.33)2 in the definition (3.7)3: 

Ma = M^ (T = 0) = [ L Pin 2 12 + U ) dV ] T= . (3.16) 

We may express the velocity of the PN mass center (3.8) as an average of the velocity inside the 
body a. Applying (3.4) with w = u ( i), and first with (p = p ( i), then with q> = p, we get by using first 
(3.11), then (2.27): 

M^a (1) = p (1) u (1) dV + 0(e 7 ) , (3.17) 

M a a = L {p U(i) + x div[p(u ( i) - u)] } dV = p u dV + .L xp(u (l) - u).n dS. (3.18) 

The surface integral is merely 0(e 5 ) by (2.15), but it is utterly negligible with our definition of 
« enlarged » bodies enclosing both the perturbed and the Newtonian position of the massive part. 
Thus we may consider that 



Therefore, we have 



M a a = L p u dV [= ord(e 3 )]. (3.19) 



M^a (1) = M a a + M a l a /c 2 , (3.20) 



l (D 

where the Newtonian term is given by (3.19), and with [using also (3.17)] 
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Ma a i/c 2 = L [(pi u + p Ul )/c 2 ]dV + 0(e 7 ). 



(3.21) 



4. Equations of motion for the mass centers: i) General form of the PN equations 

To find the equations that govern the motion of the mass centers, we simply integrate the local 
equations of motion, using the following result [deduced from Eq.(3.3)]: 

d T z i + d i (z i v j )=f : => ^-j zdV = J tdV+L z[(w-v).n]dS. (4.1) 

If (O a moves with the exact velocity: w = v = u exact , then, with z = y/u, we obtain from (2.6): 

J ffla wu l dV y \ aa { ¥ fg l -g^p^-r jk¥ u^u k - t l k¥ u k )dV , (4.2) 

in which all fields are the exact fields. Coming back to the PNA and to the notations used since Eq. 
(2.20), (d a moves with the PN velocity: w = v = U(i>, and we take z = U(i) [where V(i) = P + W\/c 2 
is the 1PN approximation of y/\. Using the PN equations of motion (2.29) and (2.30), we get just 
the PN expansion of (4.2): 

dV } = L (f + fi /c 2 )dV, (4.3) 
with 

f l = pUi, (4.4) 

fx 1 = 2 kijpj + y/i U i + pAj- 2UpU i - pu j u k - pu j d T k i} . (4.5) 

(The terms , and pi, , have been removed from/' and// respectively, by using the divergence 
theorem and the definition of « enlarged » bodies, see Sect. 3.) As is obvious in the varying units 
(see §2.3), we may identify the powers of 1/c in Eq. (4.3). Owing to (2.24), we may then insert 
Eqs. (3.19) and (3.21) respectively. We obtain thus 

M a a l =LpUidV, (4.6) 

M a l ti x l + i ai =LfidV, (4.7) 

where 

I ai = [p + p (u 2 /2 + n+U )]u 1 dV. (4.8) 
To simplify the calculation of the PN correction (4.7) to the motion, we rewrite (4.7) as 

Mati^ + i ai =J ai + K ai , (4.9) 

with 




a 

AT 



L P u i + (pu 1 i + y/ 1 u i )/, 
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K ai = L [2 kij Pj +pU ti - 2UpU t - ^ pu j u k - pu j d T kij] dV (4.10) 

and 

J ai = L a (ViU,i + pA, ddV. (4.11) 

Equation (4.9) is the equation for the PN correction to the motion of the mass center of body (a). 
Note that f l and K 6 " depend only on the first-approximation (Newtonian) fields. To evaluate J a \ we 
adapt Fock's method [5, §§ 75-77]. If we define U* = U + A/c 2 and t = <7(i) = p + oyc 2 , we have 
Ao x /c A = 0(e 8 ) by (2.15), hence 

L tU*, i dV =L a pU, i dV + J ai /c 2 + 0(e 8 ) = LpU {a \ dV + J ai /c 2 + 0(e 8 ) [=ord(e 4 )]. (4.12) 

Here and in the following, we use Fock's decomposition of any field Z, integral of some density 
vanishing outside the bodies, into « self » and « external » parts z a and Z {a) : 

Z = \edV = z a + Z ia) , Za =L a OdV, Z {a) ^ h ^ a L b OdV. (4.13) 

The second equality in (4.12) is due to the fact that, owing to Poisson's equation (2.25) and the 
boundary conditions (2.35), U is the Newtonian potential associated with p (let us denote this by U 
= N.P.[p]), hence the integral of the self-force p u a , / vanishes. Defining 

W(X, T) = \ GR p(x, T ) dV(x)/2, R = \X-x\, (4.14) 

we have AW = U. Thus (2.26) may be rewritten as 

AB = - 4nGo u B=A- d 2 W/dT 2 . 

The source of B, thus 0\, vanishes outside the bodies - contrary to the source of A. Hence it is 
natural to expect that B should be the Newtonian potential associated with 0\. Therefore we impose 
on B the same conditions as Eqs. (2.35)i_2 for U, which ensures that indeed B = N.P.[(7i]. {But 
d 2 W/dT 2 does not obey Eqs. (2.35) [see Eq. (B6) below], thus A does not behave like a Newtonian 
potential. } Thus we get: 

d 2 W r G / v d 2 W 

A(X, T) = B + — r = -o x x,r dV(x) + — Y . (4.15) 

Hence U* = U' l ' + (d 2 W/dT 2 )/c 2 with f/ f = N.P.[t]. Therefore, we have as for (4.12) 2 : 
L tU* tdV =L rU Ua) idV +L T(d 3 W/dx i dT 2 )/c 2 dV 

a a ' a 

= L a TU Ua \idV +L a p(d 3 W/dx i dT 2 )/c 2 dV+0(£ & ). (4.16) 
By the definitions of Tand Jj \ we have 

U 1 = U + B/c 2 , fl = N.P.[(7i]. (4.17) 
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From (4.12) and (4.16), it thus follows that 

J ai = Lo l U (a \ i dV +L a pB (a \idV +Lp(d 3 W/dx i dT 2 )dV+0(£*). (4.18) 



5. PN corrections for the mass centers: ii) well-separated rigidly-moving bodies 

Since the two first integrals in Eq. (4.18) involve only external gravitational fields, we may use 
« multipole expansions » : for x e (0 b and Xe (0 a with a £ b, 



1 1 



1 (x-b).(X-b) 



R X-x X-b 



X-b 



+ Gt)> r] = Sup a ^ b (r b /\a-b I) (5.1) 

K 



[r b is the radius of body (b)] . However, the source of B includes a self-field: we have 



whence 



Oi = pi + p (u72 + 17+ U ) = (p a + pU (a) , % = pi+p (u 2 /2 + n+u a ) (5.2) 
5 (a) /G = T. b * a L oJR dV=Z h * a ( L (p h /R dV + L pU {h) /R dV ) (5.3) 



Using this, it is shown in Appendix A that, for the case where each body undergoes a rigid motion 
at the Newtonian approximation: 



u^d' + n {a) ji (x j -a j ) inside (0 a , (Q % + Q {a \j = 0) 
and, as does Fock, neglecting 0(rj 2 ), but not 0(rj), with respect to 0(rj °), one has 



(5.4) 



B (a) (X) 
G 



M b [-h 2 + U^ m + qb + Ml U fn^b k + q bl + Ml(b[-b l mr-b l ) , 

h} ix -bi + ix -bp 1 

(Xeco a ); (5.5) 



L pB {a) jdV = M a 



/X=a 



(5.6) 



L (TiU (a \idV = 



l 

M a [-k 2 + U (a \ a )] + g a + Mi — — + 

2 )SX X=a 



+ U fa \t } a k + g aj + M\(a{ -a j K dx J dx j 



d 2 [J (a) 



X=a 



In Eqs.(5.5) and (5.7), f a) jk is the inertia tensor: 



(5.7) 
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£ a) jk = Lp(x j - a j )(/ - a*) dV, (5.8) 

and 

C a = (87> 1 le a )/3, Cat = (5T ai + Se ai - 17J/2, (5.9) 

where T a , e a , T ai , £ ai and r] a , are introduced by Fock [5, §74] and depend only on the rotational 
and self-potential energies and their distribution: 

T a = L pn a dv= a {a) ik a (a) jk f- a \n, n a (x) = a {a) ik a % & - a'w - d)n, (5.10) 

Tai^LpGatf-a^dV, (5.11) 

£ a = L P u a dV/2 = | space (grad u a f d V/871G = B (a \ k , (5.12) 

£ ai = L P u a (x l - a) d V/2 = | space (grad u a f (x' - a) dV/SnG = B (a \ kk , (5.13) 

%i = (B (a \ ik + B (a \ kk )/2, (5.14) 



with 
and 



B ia) u = Lpace [5«(grad u a ) 2 /2 - u a , k u aJ ]dV/4nG (5.15) 
B ia) i u = Lpace [5«(grad Ma ) 2 /2 - Ma ,* w a ,,] (x ! ' - a ! ) dV/47tG. (5.16) 
Hence, we get finally from (4.18) and (5.5)-(5.7): 



J m = GM 



|3 



+ a a^r^ + \Pai + M 1 a al)——+L a ' 
da da d a J 



d X' x=ah ^J L \X-b\ |X-b|' 

(5.17) 

where the notation <9/eta ! is used for (d/dX')x = a when no confusion may occur, and with 

cc a =M a [a 2 /2+ f/ fl) (a)] + & + B ai ^ ^\jn (a ) k d k + M a l a l , (5.18) 

L a, = L p(^V/^Jc ! '<97^)dy. (5.19) 

Thus, if the bodies are « well-separated » [r\ « 1 in Eq. (5.1)] and rigidly moving in the 
sense of Eq. (5.4), then the equation (4.9) for the PN correction to the mass-centers motion depends 
only on Newtonian quantities, and, linearly, on the PN corrections ai to the positions of the mass 
centers and the (constant) PN corrections M a l to the masses. Note that the ai 's are precisely the 
unknowns of Eq. (4.9). Clearly, the two assumptions above are well-justified in celestial mechanics, 
especially in the solar system (and the more so as we are here calculating PN corrections). The 
Newtonian quantities involved are: M a , a, a, I {a )i, Q C, a and Q , plus I ai [Eq.(4.8)], K ai 
[Eq.(4.10)], and L m . Theoretically, the Newtonian fields, hence all these quantities, should be 
known at the stage of PN correction, in which case the resolution of (4.9) would be very easy. 
However, it remains to obtain explicit expressions for l a \ K m and L a \ by using the separation into 
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« self » and « external » fields; a first step for this is done in Appendix A, leading to Eqs. (A10)- 
(All), (A13), (A21) and (A24). It appears from these equations that explicit expressions seem 
unlikely to be obtainable unless further simplifications are introduced. 

6. PN corrections for the mass centers: iii) very-well-separated spherical bodies 

6.1 Introduction of the « very good separation » and « spherical symmetry » assumptions 
In the solar system, the separation between bodies is much larger than their dimensions: e.g. the 
ratio (Earth radius)/(Earth-Sun distance) is rjE-s = 4.10 -5 , even the ratio (Sun radius)/(Earth-Sun 
distance) is still r] s . E = 5.10" 3 . Now, in Eqs. (All) and (A24) for the « external » part of the PN 
corrections, the second term on the r.h.s. is systematically of order 77 times the first term. More 
precisely, in (A 1 1), we have 

M a d 1 U (a) (a) ~Y, b ± a GM a M b al\2L-h\, (6.1) 
f a) jk Q {a) ki U {a) j(a) = Z b * a M a r a 2 x(Vj a) /r a ) x GM b /\a - bl 2 , (6.2) 
(with V ro t (a) = & {a) r a the linear rotation velocity), hence the ratio is of order 

Pat « f\ab Vj a) ld . (6.3) 

Thus, ri a h is multiplied by the other parameter 

/ a) =Vj a) /d, (6.4) 

for which estimates shall be given below. In the same way, on the r.h.s. of (A24), we have two 
sums. The two terms in the first sum are of the same order of magnitude. The first term in the 
second sum is of order (X <a) ) 2 times the former terms, and the second one is of order p ha times those 
same terms. Now recall that the equations for PN corrections, as well as the exact equations of the 
scalar theory [and in contrast with the first-approximation (Newtonian) equations], must be written 
in the preferred reference frame. From what we know about intra- and inter-galactical relative 
velocities, as also if we make the working assumption that the cosmic microwave background is at 
rest in that « ether », we expect that the velocity of the solar system through the « ether » should be 
something like 300 km/s, which would imply that the absolute velocities of the bodies would all 
have this same order of magnitude. Hence the parameters x ia \ instead of being of the order of unity 
(for Jupiter or Saturn) to 10 (e.g. for the Earth) in the solar system, as if one would take the orbital 
velocities in Copernic's reference frame, should be of the order 1/30 to 10 . Thus, the ratios p ah 
should be of the order 10" 4 to 10" 8 , i.e. very small, and (rf a) ) 2 should be of the order 10" 3 to 10" 6 , 
i.e. quite small also. Yet it seems dangerous to assume in advance the absolute velocity, hence we 
shall not neglect the self-rotation of the bodies, since for giant planets it contributes a kinetic 
energy of the same order of magnitude as does the orbital motion. 

However, we shall use the simplifying assumption of spherical symmetry for all bodies in 
the system, in order to get tractable calculations for F?". Obviously, this is not too far from the 
reality as regards the Sun, planets and satellites, but it is difficult to assess directly the numerical 
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error which is done with this assumption. What we know is the geometrical flattening at the poles, 
f a = (r a - r' a )/r a (with r a and r' a the equatorial and polar radius; the flattening is due to the self- 
rotation). This is negligible for Venus and Mercury, very small for the Earth and March (0.003 and 
0.005 respectively), but it is as large as 0.06 for Jupiter and 0. 1 for Saturn. As to the Sun, there has 
been a controversy: according to Dicke & Goldenberg [28], one would have/s « 43xl0"7959 « 
4.5xl0" 5 , whereas Hill et al. [29] found / s « 9xl0" 3 /959 « 9.4xl0" 6 [4]. The controversy arose 
because the value found by Dicke & Goldenberg would have contributed some 4" by century to 
Mercury's perihelion shift, i.e. a 10% fraction of the 43" « residual » shift. Of interest to us is the 
fact that 4" make a fraction 4/(538xl0 6 ) ~ 7.4xl0" 9 of the Newtonian main term (here the orbital 
angular displacement after one century). Note that this 4" contribution is indeed the Newtonian 
(first-approximation) correction to the orbital motion of Mercury, due to the oblateness of the Sun 
(according to its measurement by Dicke & Goldenberg). Let us provisionally admit that the effects 
of oblateness depend linearly on the oblateness f a (which is true for small oblatenesses) and that 
moreover they imply the same relative error, whether in any Newtonian term or in the 
corresponding PN correction (which is a rough conjecture). Then we expect from this example that 
even an oblateness of 0.06 (that of Jupiter) would give a relative contribution of merely some 
0.06x7. 4xl0~ 9 /(4.5xl0~ 5 ) ~ 10~ 5 to a given PN correction. Of course, this is a very rough estimate, 
yet it indicates that the assumption of spherical bodies should involve only a small error when 
calculating the PN corrections in the solar system. Anyway, this is the assumption made in GR also 
when calculating the PN corrections in the solar system [4, 18, 22, 23]. Thus, we shall assume 
spherically symmetric bodies henceforth. This assumption will concern merely the self -fields of the 
first approximation, as used for the calculations of second-approximation corrections. Therefore, 
the incompatibility between exact spherical symmetry and the presence of self-rotation is not 
relevant here. Moreover, we shall retain only the first non-zero term in the expansion with respect 
to the separation parameter J]. This is consistent with the orders of magnitude involved: if one 
would wish to account for the next order in 77, one should also remove the sphericity assumption. It 
turns out that, when one encounters two successive powers of t], the ratio of the second term to the 
first one is often of order r\tf a \ Recall that x (a) is small in Copernic's frame, except for the giant 
planets, and that rf a) should be small for all bodies in the preferred frame. We conclude that the two 
foregoing assumptions involve a relative error which is almost certainly smaller than 10 , and 
which is actually likely to be smaller than 10 , when calculating the PN corrections in the solar 
system. We emphasize that these two assumptions, together with the additional neglect of the self- 
rotation (which we shall not impose), are made when calculating PN corrections for the solar 
system in GR, since the Einstein-Infeld-Hoffmann equations (based on point masses) are very 
generally used for practical calculations [18]. Since, in GR, Copernic's frame is relevant, the 
neglect of the self-rotation of giant planets would need some justification. 

6.2 Simplifications resulting from spherical symmetry and « very good separation » 
In (All), the second term on the r.h.s., being 0(77) times the first term, now disappears. In 
(A24), the third and fourth terms disappear also, because the inertia tensor of a spherical body is 

t a \ = Ya Sij, 7a = (4%/3)\ ra r 4 p(r) dr = J (fl) r 2 p dV/3, (6.5) 

which is a constant for the case of rigid motion. Passing to Eqs. (A 10) and (A21) that are concerned 
with the « self » part of the PN corrections, we note that the ratio 
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a^ ( V a) ^ (a V(^a 2 a'), 

which is of the order (X (a) ) 2 , is like the ratio of the rotational to translational kinetic energies, hence 
the numerator term above should be conserved in (A10); with (6.5), it becomes: 

d k Q ( W2 W,. = 7a i2 a {a) jk d \ (6.6) 

In (A21), the ratio 

which is a priori (independently of any shape factor) of order x (a \ cancels in the case of spherical 
symmetry, currently considered. In fact, integrals like £ a j and T aj [Eqs. (5.13) and (5.11)] are exactly 
zero in the more general case of « axial plus planar » symmetry where the field to be barycentered 
(pu a or pQa) is invariant by rotation around some axis and by symmetry with respect to some plane 
which is perpendicular to that axis - the axis intersecting the plane at a. Note that this is very nearly 
the case for most massive celestial bodies. Moreover, for spherical symmetry, we have [using Eq. 
(C4) in the definition (5.15)]: 

B (a \ k = 8 ik \~ 4tt[(1/2) - (l/3)]r 2 (du a /dr) 2 dr/4%G = 5* £ a /3 (6.7) 

and 

B ia) i u = J"~ [5 k il2 - n k ni]r\dujdr) 2 m dr d(o/A%G = 0, (6.8) 

because 

/sphere «( dCO = 0, I sphe re «( tlj tl k dCO = 0. (6.9) 

Compiling Eqs. (AlO)-(All), (A21) and (A24), we obtain under our assumptions: 

f = [M a a 2 12 + 2T a + 4e a ] d'+jaQ % Q {a) jk d k + M a d l U (a) (a), (6.10) 



L m = -{2/i)e a d 1 -— ^Y.M b 

1 b*a 



^ - k d 2 \a-b\ . k .j ^ 3 |a-b| 

— ■ — Y~ b b — i k ■ 
dada dada da 3 



(6.11) 



To the same approximation, i.e. neglecting r\ as compared with 1, we may also rewrite (5.17) as 

J -L ai =a a ——+GM a — L^rr. (6.12) 

(Note that C, ai is entirely negligible as compared with la - bl, because their ratio is of order r] 
independently of any shape factor, while C, a i is exactly zero in the case of « axial plus planar » 
symmetry that applies very accurately to massive bodies of the solar system.) Now there remains 
the integral I?" [Eq. (A13)]. We note that, inside body (a), the /% tensor [Eq. (2.31)] may be 
approximated as 

kij = Uh (a \j [1 + 0(77 2 )], H a) ij J^i^, (6 .i3) 

u a,k u a,k 
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because, inside (a), we have 



VU=Vu a [l+0(rf)]. 



(6.14) 



Under the current assumptions, the calculation of the integral K m [eq. (A 13)] does not involve any 
difficulty: the components of this integral are given by Eqs. (C10)-(C13), (C23). Taking benefit of 
these, and of Eqs. (6.10)-(6.12), we rewrite the equation for the PN correction to the motion of the 
mass centers, Eq. (4.9), as: 



! , dU (a) 

M a a x l = a a ——+GM 
a a 



a dr 



X=a \ b * a 



X-b 



(2/3)£ a d l 



GM, 



Em, 



b±a 



pl^^y. 9 la-bl 



V 



da 1 da 



da'da k da J 



-3e a C/ a) >; (a) 



+ (d/dT) { -a \2e a + M a [/ a) (a)]/3 - [M a a 2 12 + 2T a + 4eJ d '' - 7« Q *2 % d k - M a d 1 U (a) (a) } 
+ [M a l/ a \n) + 2£ a ] Q^ji d j -(y a /5)[ri a) mj n {a) mj [/<(a) +2i2 (a) i7 t/ fl) ,/(a)]. (6.15) 
The PN correction to the mass is given by (A25), and the spherical symmetry of each body gives us 



U (a \x) = 'L b ^ a GM h /\x-b\ 



(6.16) 



Thus, Eq. (6.15) may be rewritten more explicitly as 



Ma'df = G 



dx' 



y a ab 



X=a\ b * a 



X-b 



GM 



b±a 



d\-\>\ - t . ; d\-\>\ 



da 1 da 



-b k b j 



d a l d a k d a-* 



- (d/d7)[fl ! '{ (16/3)e a + 2T a + M a [(a 2 /2) + (4/3)[/ a) (a)] }+ y a i2 £2 {a ) k a k ] 

+ [M a [/ a) (a) + 2e a ] i2 (fl) ;7 d j - {y a l5)[Q {a) mj Q {a) mj [/^(a) +2i2 (a) i7 i2%- ^(a)], (6.17) 



where 



<Xab = M a M b 



^±^ + £/ (fl) (a) + £/ w (b)' 



2 8 118 

+ ^*(^ + -e a + j7'fl)+Af a (Afi + -e fc + -r fc ).(6.l8) 



11 8 
-£ b + - 

3 3 



When the (standard) PNA of GR is used for extended bodies, rotational terms do also occur, of 
course. Those include coupled terms like (6.6), that involve both translation and rotation velocities 
[5, Eqs. 78.04 and 78.06] - but these terms are neglected in practical calculations based on the 
Einstein-Infeld-Hoffmann equations [18]. As emphasized at the end of §6.1, this neglect does not a 
priori seem to be well-justified in the case of giant planets. 

We shall not attempt here to further simplify this equation of motion, since anyway it must 
be entered in a computer program and this is clearly feasible: all quantities on the r.h.s. of (6.17) are 
known from the first- approximation calculation. Thus the PN corrections to the mass centers, ai(7), 
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shall be obtained just by (a double) time integration. However, we have to discuss the way in which 
the first-approximation quantities should be obtained. 

7. Link with observational data. Remarks on the algorithm for numerical calculations 

7.1 General remarks on the link betwen observation and prediction in a theory of gravitation 
Let us consider some « relativistic » theory of gravitation, (T), for which one may define a precise 
Newtonian limit in terms of a parameter equivalent to the parameter X = e 2 = (1 -/ m in)/2 in the 
scalar ether-theory [see Eq. (2.11); we have seen since that indeed the square e 2 is more relevant, 
because the scalar theory admits expansions in powers of e 2 ]. More precisely, let us assume that the 
equations and fields of theory (T) admit first-order expansions in terms of X, with the zero-order 
equations reducing to NG, and with the remainder terms 0(A 2 ) being negligible as compared with 
the observational error (note that this is currently the case as regards astrodynamics in the solar 
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system, with £q ~ 10 ). Let D = (Dj) j = ^ . , j be a representative set of observational data 
(Keplerian parameters of planetary orbits and their time-drift, etc.). In order to perform the 
theoretical « prediction » of these data, theory (T) introduces some finite set of numerical 
parameters, a. Let a = (a, ); = i, / o be the subset which is needed for the zero-order calculation 
(this will include the first- approximation masses M a and some higher-order multipoles of the first- 
approximation mass density, etc.). Let a 1 = (a, 1 ), = i, ^ be the complementary set which may be 

needed, in addition to a (or rather to a part of a ), for the first-order (second-approximation) 
correction. Thus 



and the theoretical prediction for data Dj is the sum of the zero-order (« Newtonian ») calculation 
and the first-order (PN) correction: 



(Of course, the functional expressions Fj are not analytical ones and in practice will involve 
numerical calculations.) 

Now it should be realized that the set of parameters a can hardly be accessed independently of the 
set of observational data D. For instance, it is clear that one cannot weigh the planets and that the 
masses M a are free parameters of the model [i.e., of the first-order approximation of theory (T)]. Of 
course, some parameters such as the velocity of light (which does enter the expressions for PN 
corrections) can be accurately measured independently of any celestial mechanics, and we just 
remove them from the list a. Thus, the process of theoretical « prediction » involves essentially a 
fitting of the observational data Dj by the functional expressions Fj. In practice, this will be a least- 
squares fitting, so one will have to search for the values of the free parameters a that minimize the 
residual: 



a = (a , a 1 ) 



(7.1) 



Dj theory = Fj (a) = Fj (a ) + XFj 1 (a). 



(7.2) 



R(a) = Zj [Fj (a) - Dj] 2 = Min. 



(7.3) 
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But, in just the same way, what is called « Newtonian celestial mechanics » is essentially a fitting of 
the data Dj by the first- approximation expressions Fj °, and this gives the Newtonian values a N = 
(ccP)i = i, / of the zero-order parameters: 

R (a N ) = Y,j [Fj (a N ) - Dj] 2 = Min. (7.4) 

The Newtonian values of the zero-order parameters have no reason to be optimal for the PN 
calculation, simply because the functional expressions to be minimized are different. In particular, 
the Newtonian masses M a N have no reason to give adequate values of the zero-order masses M a for 
PN calculations. 

To get some idea about the difference that makes, it is reasonable to assume: (i) that « the 
theory (T) is good » in the sense that its PN approximation (with the optimal set of parameters a) 
gives predictions for D which coincide with D up to an 0(X 2 ) error [thus, the improvement of NG 
by theory (T) should not be found merely in the solar system, but instead should be verified each 
time one tests a weakly gravitating system] ; (ii) that « NG is not bad » in the sense that the 
difference between the Newtonian values a, N of the zero-order parameters (values that are optimal 
for the first approximation) and their correct values a? (that are optimal for the second 
approximation) is small; and (iii) that, however, « the theory (T) is necessary » in the sense that the 
difference between the Newtonian prediction F,°(a N ) and the observed data Dj is ord(A). {We adopt 
the varying units of mass and time, [M] e = ^[M] and [T] e = £~ l [T] (§2.3), so that D, a, F°(a°) 
and F\a) are ord(A ).} From (i) and (ii), we get by (7.2): 

Fj (cc N ) + (dFj °/da?) 5a? + XFj 1 (a) = Dj + 0(A 2 ), da? = a? - a, N . (7.5) 

Hence, by (iii), the quantities 

(dFj °/daP) da? = [Dj - Fj (a N )] - XFj 1 (a) + 0(A 2 ) (7.6) 

are 0(A) anyway, and should a priori be ord(A), unless miraculously the a N were systematically 
optimal zero-order parameters for the first-order calculation. Since dFj °/da^ = ord(A°) in the 
varying units, this means that, in the same way, da? should be ord(A) [and is 0(A) anyway]. If, in 
the PN fitting (7.3), one assigns the Newtonian optimal values 0C N to the zero-order parameters, one 
shall get values a' 1 ^ a 1 and so the whole set of PN parameters shall be a' = (a N , a' 1 ) ^ a. Then, the 
total PN prediction obtained so will differ from the optimal PN prediction by: 

8Dj = Fj (oc N ) + XFj 1 (a') - [Fj (a ) + XFj 1 (a)] = (dFj °lda?) da? + X[Fj 1 (a') - Fj 1 (a)] (7.7) 

and so, a priori, will be ord(A). Thus, the error made in taking the Newtonian optimal values 0^ for 
the zero-order parameters is a priori of the same order [ord(A)] as the PN correction. Admittedly, 
we cannot exclude the case that a « miraculous » cancellation would make the difference (7.7) 
0(A 2 ), hence acceptable for the first PNA, but in the absence of any proof of that, we cannot assign 
the Newtonian optimal values a N to the zero-order parameters a . This applies to any theory of 
gravity admitting asymptotic expansions for weak fields and whose zero-order approximation 
reduces to NG. 
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7.2 Relevant parameters in the scalar theory and the question of the preferred frame 
Let us list the parameters needed to perform the calculation of the PN corrections as given by Eq. 
(6.17). These are the first-approximation masses M a , and other integrals of the first- approximation 
fields: T a [Eq. (5.10)], £ a [Eq. (5.12)], y a [Eq. (6.5)] and Mj [Eq. (A25)]. Note that the presence of 
these integrals means unambiguously that the internal structure of the bodies influences the motion 
at the level of the PN corrections. Moreover, in (6.17), we have functions of the absolute time T 
(i.e. the inertial time in the preferred reference frame E), namely the first- approximation positions, 
velocities and accelerations of the bodies, a(7), a (7) and a (7) and the rotation velocity tensors £l (a) 
for a = 1, TV (actually, the latter tensors may be considered constant for the purpose of PN 
corrections at the accuracy envisaged here). But the accelerations are given by Eq. (4.6), hence are 
determined by the masses M a , the inertia tensors I (a \ [and if desired some higher-order multipoles 
Ihi g her (a) of the zero-order mass density p a ,] plus the current positions a(7). To determine the zero- 
order problem, we also need the initial velocities and positions. Thus, the list of the zero-order 
parameters is as follows: 

a = (M a , I (a) , [Ihi g her (a) o,] a , a o, (a = 1, AO), (7.8) 

where index means initial value. (Since the bodies are assumed to have a rigid motion, the inertia 
tensor I (a) [and Ihi g her (a) ] merely rotates with the rotation velocity & a) .) 

The positions, velocities and accelerations have to be given in the frame E, yet the astronomical 
observations, e.g. in the solar system, are expressed in some moving frame F - such as Copernic's 
reference frame C, whose axes are defined with respect to the « fixed stars », and whose origin is at 
the global mass center of the solar system. In order to properly define the correspondence between 
observables in the different frames involved, say E and F, we shall assume that F undergoes a pure 
translation, at some velocity V(7), with respect to the preferred frame E. This means that, as seen 
from E (and, in particular, using the « absolute simultaneity », defined by the Poincare-Einstein 
synchronization applied in the frame E), each observer in F has the same velocity V(7). For 
observations in the solar system, the relevant absolute velocity V shall be that of its mass center, i.e. 
one will have F = C. For PN corrections in the solar system, the relevant time scale is the century 
[for larger time intervals, one should « reinitialize » the problem, see around Eq. (3.14)]. Therefore, 
we shall assume here that V is a constant, so that the correspondence between E and F is a Lorentz 
transformation of the flat metric 7° [14] - a special one, if one takes the Ox 1 and O'x' 1 axes parallel 
to V: 

x' 1 = Yv(x l -VT), T = y v (T-Vx 1 /c 2 ), x' 2 = x 2 , x' 3 = x 3 . (7.9) 

Hence, relevant physical fields in the moving frame (now denoted E v ) are the transported proper 
fields: 

P'exact (T\ X')= Pmct [T(T , x'), X(T\ x')] (7.10) 

and the like for p* and 77, plus the Lorentz-transformed velocity and acceleration: u' exa ct = dx'/dT" 
and A'exact = du' exa ct/dr'. As to the gravitational field, we shall consider the transported field 

f'(T',x')=f[T(T',x'),x(T',x')], (7.11) 
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though it has not a direct physical significance. Since the fields in the frame E admit asymptotic 
expansions, we may obtain ones for the fields in the frame E v . Because we must assume that V/c = 
0(e) [14], it follows that the Lorentz transformation (7.9) involves the small parameter. Therefore, 
we have to account for the expansion of (7.9). Neglecting 0(e 4 ), the latter is: 

x' 1 = [1 + V 2 /(2c 2 )](x l - VT) [1 + 0(e 4 )], T = {T+ (V 2 /c 2 )[(T/2) - (x/V)] }[1 + 0(e 4 )], (7.12) 

where it has been used the fact that x is in the near zone, i.e. x = 0(e ). In particular, when one 
neglects 0(e 2 ), it becomes the Galileo transformation: 

x' = (x- Vr)[l + 0(e 2 )], T = T[1 + 0(e 2 )]. (7.13) 

We now have to insert the initial expansions into the expressions of the transported fields, e.g. 
insert (2.17) into (7.10), and to account for (7.12) and (7.13). This provides an expansion for the 
transported fields: e.g. 

p'aUT, x')=p ex , ct [T(T\ x'), x(F, x')] = [p(T, x) + Pl (T, x)/c 2 ][l + 0(e 4 )] 

= [p\T, x')+p\ (T\ x')/c 2 ][l + 0(e 4 )]. (7.14) 

Due to (7.13), the zero-order term of the transported field is simply the Galileo transform of the 
initial zero-order term, i.e. 

p\T',x')=p(T',x' + \T) (7.15) 

so that the notation p' is unambiguous. The same is true for p* and 77, of course. In the same way, 
we have U'(T\ x') = U(T\ x' + \T) since, from (7.11) and (7.13): 

/' (T\ x') =f(T, x) = 1 - 2U(T, x)/c 2 + 0(e 4 ) = 1 - 2U{T\ x + YT')/c 2 + 0(e 4 ). (7.16) 

Moreover, it is easy to show [14] that the Lorentz-transformed components y'/jy of the physical 
metric y satisfy 

r'oo(7", x') = 1 - 2 U\T\ x')/c 2 + 0(e 4 ), (7.17) 

so that the transported potential U ' has a direct physical significance. As to the velocity and 
acceleration, they were not defined by a mere transport, but we get easily for the zero-order terms: 

u' (r, x') = u(r, x' + vr) - v, A' (t\ x') = du'/dr = (du/dr) (r, x + yt). (7.18) 

The equations for the exact fields in the moving frame Ey are obtained by substituting for the fields 
in the preferred frame E their expressions in terms of the fields in Ey, like (7.10) and (7.11), in Eqs. 
(2.3), (2.6) and (2.7). Using the Lorentz transform, one expresses the derivatives in E as derivatives 
in Ey [14, Eq (49)]. As one may easily check, it follows from these definitions that the zero-order 
fields in the frame Ey obey the equations of NG, Eqs. (2.25), (2.27) and (2.29) with primes. This 
was expected since (i) the zero-order fields in the frame E obey NG, (ii) the zero-order fields in E 
and in Ey exchange by Galileo transformation, and (iii) NG is Galileo-invariant. 



22 



7.3 The gross structure of an algorithm for celestial mechanics in the scalar theory 
Having found a firm basis for NG in the scalar theory, we may use it, in particular, to assess the 
integrals T a , e a , y a and M a l [Eqs. (5.10), (5.12), (6.5) and (A25)]. Thus, in the solar system, we may 
use Newtonian estimates of the density profile p = (p a (r) for the different bodies (a) [30-31] and 
compute these integrals. As emphasized in §7.1, such purely Newtonian calculations of the first- 
approximation quantities are correct only up to an 0(e 2 ) error, but this does little matter here, 
because £ a , E, a and a are used only in second-approximation (PN) corrections, hence an O(^) error 
on them makes an 0(e 4 ) (third-approximation) error in the total PN estimate. Thus, the PN 
parameters a 1 reduce to the absolute velocity of the moving frame Ey used for observations, i.e. to 
the vector V. In practice, that moving frame could be Copernic's frame C. However, most 
observations are performed, not in a such uniformly moving frame directly, but instead in local 
frames with a more complex motion, e.g. with Earth-based telescopes. In GR, formalisms exist to 
pass from « local » frames to « global » ones [32, 33]. We shall not discuss here the way this 
correspondence has to be modified in the scalar theory and will assume that the observables D are 
expressed in some uniformly moving (« inertial ») frame Ey - but V is not a priori known. 

The skeleton of an algorithm may then be proposed for PN calculations in the scalar theory: 

0) Initialization: Perform the purely Newtonian calculation, i.e. solve (7.4) by iterating the (non- 
linear) zero-order calculation alone, and thus obtain the Newtonian values Ct N of the zero-order 
parameters (7.8). Of course, it will be more convenient to make that step in Ey. Since the list (7.8) 
is given in E, one will have to Galileo-transform the zero-order parameters from Ey to E. As it 
appears from (7.8), only the velocities ao shall be affected: they shall be modified according to 
(7.18)l 

1) a) Some estimate for the zero-order parameters a and for V being known, compute the zero- 
order motion (as was done for each iteration of step 0), and compute the PN corrections to the 
motion by (6.17). Then Lorentz-transform the total motion, to 0(l/c ) accuracy, back to Ey. 
Estimate the residual (7.3). 

b) Using a minimization algorithm, make loops on a) to optimize a and V. 

Two different starting estimates seem possible for V: either 0, or the velocity against the 
cosmic microwave background. The decisive points for the theory shall be (i) the accuracy, i.e. the 
level of the final value of the residual (7.3); and (ii) the magnitude of the optimal velocity vector V. 
Since the theory predicts exact Schwarzschild motion in the spherical static case, one is tempted to 
conjecture that a correct accuracy (improving significantly NG) should be reached. But, for which 
vector V? If V = IVI were found too small (some km/s), that would mean Copernic's reference 
frame is nearly at rest in the ether - which seems difficult to justify. On the other hand, if a good 
agreement were found for a significant value of V, that would be a strong argument for an ether. 
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Appendix A. Integrals I m , J m , K ai and L m and Ma (well-separated rigid bodies) 

To evaluate integrals which are similar to I m and J m , Fock [5, §§ 76-77] uses the relations 

pn-pu a +p = pn a , (Al) 

u 2 /2 = a 2 12 + d k n (a) jk (^ - d) + Q a , (A2) 

valid inside body (a). Equation (A2) follows from the particular form (5.4) assumed for the velocity 
field, and Eq. (Al) is deduced from the Newtonian equation of motion (2.29), accounting for (5.4) 
and for the expression of the elastic energy in a barotropic fluid: 

p*{q) p*(p) 
Let us evaluate J ai . It follows from (Al) and (A2) that 

p (u 2 /2 + 77 + u a ) = p [a 2 12 + d h & - d) + 2Q a + 2 u a ] - p. (A4) 
Using the definitions (5.10)-(5.14) and Fock's integrals (74.24) and (74.25): 

3L pdV=£ a -2T a , (A5) 
2Lp(x i -a)dV=ri ai -T ai , (A6) 
we get then [remembering the definition (3.9) of the mass center] 

L P (u 2 /2 + 77+ u a ) dV = M a a 2 I2 + (ST a + 1 le a )/3 (A7) 

and 

L P (u 2 /2 + 77 + u a ){x [ - a) dV=d k G ( V^; + (5^ + %e ai - 7] ai y2. (A8) 
We have from (5.1): 

L b (pAX-x\ dV(x) = [L b (p dWIX - bl + {X 1 - b) L b (faW - b')dV(x) AX - bl 3 ][l+ 0(J] 2 )]. (A9) 

We apply this to <p = (p h + pU (b) and use (5.3). For (p h , we use (A7)-(A8) and (3.9) 2 ; for pU (b \ we 
make a Taylor expansion of lf b) at b and use (3.9) i. This gives us Eq. (5.5). Making a Taylor 
expansion of B ia) at a and using (3.9)i, we get (5.6). Equation (5.7) is obtained in the same way. 

To calculate f [Eq. (4.8)], we substitute the expression (5.4) for u in (4.8) and use (A5)- 
(A8) to obtain 

L [p + P (u 2 /2 + 77+ u a )] u l dV = [M a a 2 12 + 2T a + 4e a ] d' + id ^ + 2T aj + 4e aj )Q (a) y7 .(A10) 

Making a Taylor expansion of U ia) at a and using (3.9)i, we get 
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L pU {a) u [ dV = [Ma a ' f/ (a) (a) + 7% Q ia) ki t/ (fl) ,/a)] [1+ 0(r? 2 )], (Al 1) 

and I m is obtained by summing (A 10) and (All). 

As to K"': using (2.3 l)i and applying (4.1) with z = pKjd , neglecting \da> a Z Ui.n dS/c 2 (as is 
plainly justified, see Section 3), we find by (2.29) that 

L dr l j k pu j u h + pu j d T kij) dV = L pu j [(d T kij + k u u k ) + (k ikJ - k jk j) u k ]dV 

= (d/dT) (L p kiju j dV) - L ky (pUj - Pj ) dV + L pu j (k ikJ - k jkJ ) u k dV. (A12) 

Substituting this in the definition of K m (4.10), we get: 

K ai = L a kijipj + pUJdV + L a pU AV + L a (~2UpU ddV- (d/d7)(L a pk ijU j dV) + L a pu j u \k jKr k ikj )dV 
Ki d + K 2 ai + K 3 ai - (d/dT)K' ai + K A ai . (A 13) 

We decompose L m into « self » and « external » parts: 

L ai = L p (d^a/dx 1 d T 2 ) dV+ L p (d^Vdx'dT 2 ) dV. (A14) 

The transformation of the « self » part is rather involved. From (4.16), it follows that 

I 

„ tu* a , i dV = LtU "'" (a \ i dV + [L p (d 3 \^ a) /d x i dT 2 )dV + L P(d 3 w a /d x' d T 2 ) dV ] /c 2 

-L TU* {a \idV+0(e\ (A15) 

and since U* (a) = U W) + (dV a V dT 2 )/c 2 , we get 

L p (d^Wa/dx'dT 2 ) dV/c 2 = L a ni* a , i dV + 0(£ 8 ). (A 16) 

It is proved in Appendix B that 

L TU* a jdV= I space U* a , ,■ (d 2 U a /dT 2 ) dV ] /(4%Gc\ (All) 
Ispace U* a , (d 2 U a /dT 2 ) dV = I space U a , (d^/df 2 ) dV + 0(£ 8 ), (A18) 

Lpace u a , i (d^JdT 2 ) dV = (d/dT) L 4;iGp (d 2 w a /dx l dT) dV . (A19) 

From (A16)-(A19), we obtain 

L p (d^Wa/dx'dT 2 ) dV = (d/dT) L p (d^Jd^dT) dV + 0(£ 8 ), (A20) 

which we were not able to prove in a simpler way. The integral on the r.h.s. is calculated by Fock 
[5, Eq. (76.11)], which gives 
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L P (d 3 w a /dx'dT 2 ) dV = (d/d7)[- £ a d [ + B (a> ik d k -n (a) ji£ a j + Q ia) jk B (a> jik ]. 
As to the « external » part of L ai , we get by a Taylor expansion and (3.9)i: 

L P (d^Vdx'dT 2 ) dV = [M a (^^x'^Ka)] [1+ 0(r? 2 )]. 
Using Fock's formula (76.34): 

d 2 W (a) l r -,d 2 \x-b\ lr„.(j) <? 3 |x-b| 



dx'dT ^ b^a dx l dx 4 i _ tfl dx l dx J dx 

2 

we obtain thus [neglecting 0(7] ) as compared with 0(r] ), as also does Fock to get (A23)]: 



r <9 3 W (fl) ..^ 2 |a-b| <9 3 |a-b 

b*a 



t ,. , u ~ u v - ~ II / , Mu 



da l da k d a l d a k d a 7 y 



_ 1_ y <9 3 |a-b| _ . (b) ^ <9 4 |a-b| "jl 

2^ a \ Jk da i da J da k jk da l da j da k da l P 



According to (A 14), L m is obtained by summing (A21) and (A24). 

Finally, we get by inserting (A2) into (3.16) and with the usual Taylor expansion of 
Ma = [{M a [a 2 /2 + t/ fl) (a)] + T a } T = + 2e a ] [1+ 0(r? 2 )]. 
Note that £ a is a constant as far as we assume that (a) has a rigid motion. 

Appendix B. Proof of the equalities (A17), (A18) and (A19) 

Proof of Eq. (A17). Since U* = U + A/c 2 , we have by (2.25) and (2.26): 

- 4tcGt= AU* - d^/d^/c 2 , 
whence [remembering the definition (4.13) for the self fields] 

4llG L T U* a jdV= Jspace U* a , / (d 2 U a /dT 2 ) dV/c - I S p ace W* a , / All* a dV. 

In this work, we define the integral over the whole space as the limit (when it exists): 

Ispace 9 dV = lim r _> oo I|x-AI < r 9 dV, 
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which depends, in general, on the choice of the origin A. For the integrals considered, we shall take 
A = a, i.e. the first- approximation mass center of body (a), but actually these integrals do not 
depend on the origin. We obtain thus by the divergence theorem: 

IspaceW*a ;i Aw* a dV= lim L [rijU* a<i u*aj-*iu*ajP*aji/Z\dS = lim Is r(g*.n)g*-g* 2 n/2]dS(B4) 

with g* = Vu* a and where S r is the sphere Ix - al = r. Since u* a = d 2 wJcfT 2 /c 2 + N.P.[t„], we have 
g* = V(d 2 w a /df 2 /c 2 ) + 0(l/r 2 ). Now by expanding Ix - yl in a Taylor series at x - a, we get: 

w a (x, T) = L G\x - ylp(y,7)dV(y)/2 = GMJx - al/2 + GP^Sg- n^)/(4lx - al) + (B5) 

d T w a = -GM a n.k/2 + 0(l/r), d T 2 w a = -GM a ii.a/2 + 0(l/r), (B6) 

V(d 2 w c /dT 2 ) = -GM a z(n)/(2r) + 0(l/r 2 ). [ z(n) = a-(n.a)n, n = (x-a)/r, r=lx-al. ] (B7) 

Since z.n = and since z(n) = z(-n), it follows that the limit in (B4) is zero. This proves (A17). 

Proof of Eq. (A18). As for d T 2 w a , one finds that 

u a (x, T) = L Gp(y,7)dV(y)/lx - yl = GM a /lx - al + G/ aJ !7 (3n ( n y -^)/(2lx - al 3 ) + (B8) 

d T u a = GM a ii.a/r 2 + 0(l/r 3 ), d T 2 u a = GM a ii.a/r 2 + 0(l/r 3 ) (r^oo). (B9) 

At large r, u a j is ord(l/r ) and, as we have seen, u* a j is ord(l/r). Moreover, both integrands in 
(A18) are bounded. Thus, the integrand on the r.h.s. of (A18) is (Lebesgue-)integrable, being 
0(l/r 4 ) at large r, so that the integral on the r.h.s. a fortiori exists in the sense of (B3). The 
integrand on the l.h.s. is ord(l/r 3 ), hence not integrable. However, the main term in that integrand is 

Const x [w - (w.n)n] w.n/r 3 , [ w = a ] (BIO) 

hence changes sign with n and so makes a zero net contribution in each ball Ix-al < r, thus 
admitting zero as the limit. Since the remainder is 0(l/r 4 ), therefore integrable, it follows that the 
integral on the l.h.s of (A18) also exists in the sense of (B3). Hence, the difference between the two 
terms: 

Ispace (u* a j-U aJ )(d 2 U a /dT 2 ) dV = Ispace A a j(d 2 U a /d f 2 ) &V/c = D ai (Bll) 

also exists in the sense of (B3). In the varying units introduced in §2.3, D a i is obviously 0(£ ), and 
since D ai has dimension L 6 T~ 6 , it is therefore 0(e 8 ) in the fixed units. This proves (A18). 

Proof of Eq. (A19). We have by the divergence theorem: 

Ispace u aJ d T 2 u a dV= | space [d j(d ill a d T U a ) ~ d T Uad T diUa\ dV 
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= (d/dT) LpacediUadjUadV- lim Is [d T (u a 2 )/2] Hi dS . (B12) 

We used an inversion between integration and time differentiation. In order to justify this by 
Lebesgue's theorem, we must show that: i) for any T, diUadfu a is integrable, and ii) <9r(<9,w a <9 r w a ) is 
bounded, independently of T, by an integrable function 0(x). Condition i) is satisfied, for diU a d T u a is 
bounded, and 0(l/r 4 ) at large r. As regards condition ii), we note that, due to the assumed quasi- 
periodic character of the motion (in some interval of time I), d T (diU a d T u a ) is bounded, inde- 
pendently of x and T, by some number a. Moreover, for any T, it has the form [cf. (B8) and (B9)] 

d j(d iu a d T u a )(x,7) = -G 2 Mamrij d j /\x - al 4 + R(T, x) (B 1 3) 

with 

R(T,x) = 0(l/r 5 ) asr = lx-al^oo. (B14) 

Again due to the assumed quasi-periodic character of the motion, the behaviour of the remainder 
term is valid independently of T; i.e., there are some numbers M and ro such that 

V7el R(T, x) < MAx - al 5 if lx-al>r . (B15) 

Hence di(diU a d T u a )(x,7) is bounded, for Ix - al > r , by the function 

(/>i(x) = MAx - al 5 + G 2 M a 2 ln ; n ; l{Supr e i [\d j (T)\] }l\x - al 4 (B16) 

which is integrable on the domain Ix - al > r . Thus, defining 0,- (x) = a for Ix - al < r , one satisfies 
condition (ii), so the inversion between integration and differentiation in (B12) is licit. Now, by 
(B8) and (B9), one has u a d T u a = 0(l/r 3 ) at large r, hence the limit on the r.h.s. of (B12) is zero. 
Thus: 

dV = (d/dT) I space diu a dju a dV . (B 17) 

Now we shall shew that 

iU a drUa dV . (B18) 

Together with (B17), this will prove (A19). We substitute -Au a for <\nGp and use Green's formula: 
L AnGp(d 1 w a /dx i dT)dV= lim J B -Au a (d 2 w a /dx dT) dV = lim J B -u a A(d 2 w a /dx i dT)dV + 

+ lim Is [u a (d/dn)(d 2 w a /d x ' dT)-(d 2 w a /d x 1 dT)dua/dn]dS = limi + lim 2 (B19) 

r— >°o r 

(B r is the ball Ix - al < r and S r the sphere Ix - al = r). We use the divergence theorem: 

limi = lim I B -u a (d 2 u a /d x dT) dV = lim [I s -u a (d T u a ) tii dS + I B u a jd T u a dV], (B20) 
r— >oo r r^oo r r 
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and by (B8) and (B9) we get thus: 

linii = lim J B u a>i d T u a dV = J space diU a d T u a dV . (B21) 

r— ¥°° ' 

From (B5), we find that d 2 w a /d x' d T = 0(l/r) and d 3 w a /d x l d Td x> ' = 0(l/r 2 ). Since dujd x l = 
0(l/r ) and u a = 0(l/r), we get lini2 = and, with (B21), obtain (B18). This completes the proof. 

Appendix C. Calculation of K° l (very- well-separated spherical rigid bodies) 

Let us calculate K m [Eq. (A13)] under the two assumptions made at the end of §6.1. We begin with 
K\ \ using Eq. (6.13). Since the self-potential u a of body (a) is assumed spherical, we get inside (a): 

k^x) = Uh {a \j [1 + 0(77 2 )] = [u a (r) + f/ fl) (x)] m nj[l+ 0(r] 2 )], (CI) 

where n = («/) = (x - a)/r, r = I x - a I. Since the pressure p is also spherical, we have 

I m U a Hi Hj jC, / dV = Jo U a Hi flj U: p' dV = Jo (u a p')(r) tlj dV = 0, (C2) 

because n, has zero integral on each sphere S r (henceforth, the prime will denote derivative with 
respect to r, except for K u "). In the same way, we get 

u a rii nj p u aJ dV = 0, L n pu a , i dV = 0. (C3) 

To evaluate the contribution of the external potential L^ a) (x), we make a first-order Taylor 
expansion of it at x = a and use the following formula [5, Eq. (90.19)]: 

\m nj dm = (4ti/3) dij (do = sin0 d0 d(p). (C4) 

We find thus: 

Kr = [(2/3)£ a -$a- a ] f/^, ,<a), (C5) 

where £ a is given by (5.12), and with 

$ a = - (1/3) L P (x - a).Vu a dV=- (47C/3) L P u' a r 3 dr, (C6) 
a = - (1/3) L (x - a).Vp dV=- (47C/3) L p'a r 3 dr. (C7) 
But, since Fock's Eqs. (73.16) and (73.19) [5] may be immediately rewritten as 

pu a , i + Q {a) ik Q % (x j -a i )=p J , (C8) 
we get by (5.10), independently of the sphericity assumption: 

-3^ a + 2T a = -39 a . (C9) 
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Therefore, (C5) may be written as 

K x ai = [(2/3)(e a + T a ) - 2Q lf a \ ,<a). (CIO) 
In the same way, using (C3)2 and Taylor-expanding p at a, we get by (A5): 

K 2 ai = L plf a \ i dV = (l/3)(e a - 2T a ) ,<a). (CI 1) 

We also get (to the lowest order in 77) 

K 3 ai = -2[L p u a lf a \ i dV + L P u a> i lf a) dV\ = (2% a - 4£a) ^ , ,<a) (C12) 
and, using (5.4), (6.13) and (C4), 

K ' ai = L P (u a + f/ a) ) ta rij (d j + Q (a) jk rn k )dV = d 1 [2e a + M a t/ a) (a)]/3 . (C 1 3) 

As to K 4 ai , we set y = x - a(7) and we may write using (CI): 
Kf = L a pu j u k [(u a y j y Vr 2 ), - (u a y l y k /r 2 \j]dV + L a P u j u k t/ fl) [( y j y k /r 2 )j -(y i y k /r 2 )j] dV + 

+ L a p u V[(f/ a) , ,■ y j y k - lf a \j y'y k ) /r 2 ]dV. (C14) 

Here, u j is given by (5.4), and u a = u a (r). Moreover, lf a) may be replaced by its first-order Taylor 
expansion at a, since we retain merely the lowest order in 77 ; and we have 

(J j y k /r\ i - (y 'y k /r 2 Xj = (y j 8 ik - y >8 jk )/r 2 = (nj5 ik - n^/r. (C15) 

Therefore, we still may write: 

K A ai = d j d h (I rp(r) u a (r)dr) (I (njd ik - n^dco ) 

+ (I r 2 p(r)u a (r)dr ) (f (nj8 ik - n^Q^n n,d k + Q {a \ k m d j )d(0) 

+ (I r 3 p(r)u a (r)dr ) (I (rij5 ik - H^)(i2 (a) )j m Q (a) mk n m )d(0) 

+ f/ a) (a) II rp(r) (njS* - n^) (d j + rQ (a) mj n m ) (d k + rQ {a) pk n p ) dr dco 

+ lf a \i(a) II r 2 p(r) m (nj8 ik - m8 jk ) (d j + rQ {a) mj n m ) (d k + rQ {a) pk n p ) dr do) 

+ II r 2 p(r) i(si)njn k - lf a \ j(a)ntn k ] (d j + rQ {a) mj n m ) (d k + rQ {a) pk n p ) dr dco 

= I 1 +I 2 + I 3 + I 4 + l5 + l6 (C16) 

(the r-integrals are from r = to r = r a ). With the help of Eqs. (6.9), (C4), and [5, Eq.(90.20)] 
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\m nj n k m dco = (471/15) (dijdu +8 ik dji +8 a d jk ), (C17) 
one finds in a straightforward manner that: 

Ii = I 3 = 0, (C18) 

h = (4tc/3)(I r 2 p(r)u a (r)dr)[(n {a \j d «'+ Q {a) u d j )8j, - 2Q ia) lj d j 8 u ] = 2e a £2% d j = 2e a (G) a Aa ),-, (C19) 

I4 = t/ a) (a)[| r 2 p(r)u a (r)dr](4K/3)(ri a) ji d j +Q {a )j d '-i2 (a) , 7 d J -i2 ( % d j ) = M a f/ a) (a)((O a Aa ),-, (C20) 

1 5 = (a) [(M a /3)(d l d '- <5 ;7 d j d j ) - (% l5){Q (a) mj Q ia) mj 8 il + 2Q (a) tj i2%)], (C21) 

1 6 = (M a /3)[l/ a \i (a) a V - t/ fl) , j (a)d j d ')] + (% /5)[0xf/ a) , ,<a) - OxC/^/a)] (C22) 
[y a is given by Eq. (6.5)]. Summing this, one gets finally 

K? = Walf a \a) + 2e a ] d j - (y a l5)[Q {a) mj Q ia) mj f/ fl) ,/(a) +2Q {a \j d a \ t/ fl) ,/(a)], (C23) 
or (using the formula of the double vector product) 

K 4 a = [M a l/ a \a) + 2e a ] ® a Aa + (y a /5) { 2[co a .Vf/ a) (a)] ® a - 4co a 2 Vt/ fl) (a) } . (C24) 
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